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Abstract 

Sex-biased genes are thought to drive phenotypic differences between males and females. The recent availability of high-throughput 
gene expression data for many related species has led to a burst of investigations into the genomic and evolutionary properties 
of sex-biased genes. In Drosophila, a number of studies have found that X chromosomes are deficient in male-biased genes 
(demasculinized) and enriched for female-biased genes (feminized) and that male-biased genes evolve faster than female-biased 
genes. However, studies have yielded vastly different conclusions regarding the numbers of sex-biased genes and forces shaping their 
evolution. Here, we use RNA-seq data from multiple tissues of Drosophila melanogaster and D. pseudoobscura, a species with a 
recently evolved X chromosome, to explore the evolution of sex-biased genes in Drosophila. First, we compare several independent 
metrics for classifying sex-biased genes and find that the overlap of genes identified by different metrics is small, particularly for 
female-biased genes. Second, we investigate genome-wide expression patterns and uncover evidence of demasculinization and 
feminization of both ancestral and new X chromosomes, demonstrating that gene content on sex chromosomes evolves rapidly. 
Third, we examine the evolutionary rates of sex-biased genes and show that male-biased genes evolve much faster than 
female-biased genes, which evolve at similar rates to unbiased genes. Analysis of gene expression among tissues reveals that 
this trend may be partially due to pleiotropic effects of female-biased genes, which limits their evolutionary potential. Thus, 
our findings illustrate the importance of accurately identifying sex-biased genes and provide insight into their evolutionary dynamics 
in Drosophila. 



Introduction 

Sexual dimorphism, which comprises the morphological and 
behavioral characteristics that differentiate males and females, 
is prevalent in the animal kingdom. However, despite often 
dramatic differences between sexes at the phenotypic level, 
males and females share nearly identical genomes. Thus, 
sexual dimorphism is thought to result primarily from differ- 
ential expression of genes that are present in both sexes, a 
phenomenon referred to as sex-biased gene expression 
(Connallon and Knowles 2005; Rinn and Snyder 2005; 
Ellegren and Parsch 2007). 

In Drosophila, sex-biased genes are nonrandomly distribu- 
ted across the genome (Parisi et al. 2003; Sturgill et al. 2007). 
In particular, Drosophila X chromosomes are deficient in 
male-biased genes (demasculinized) and possibly enriched 
for female-biased genes (feminized). Three mechanisms 
have been proposed for this observation. First, inhibition 
of the expression of X-linked genes during spermatogenesis 
by meiotic sex chromosome inactivation (MSG!) may result 
in selection against testis-biased genes on X chromosomes 
(Lifschytz and Lindsley 1972; Hense et al. 2007; Vibranovski 



et al. 2009). Second, dosage compensation, which in 
Drosophila is achieved by overexpression of the male X chro- 
mosome, may make it less favorable for male-biased genes to 
reside on X chromosomes (Vicoso and Charlesworth 2009; 
Bachtrog et al. 2010). Third, because X chromosomes are 
preferentially found in females, sexually antagonistic selection 
may prevent male-biased genes, which are likely beneficial 
to males and detrimental to females, from accumulating on 
X chromosomes (Rice 1984; Ellegren and Parsch 2007; 
Gurbich and Bachtrog 2008). However, despite extensive 
research in this area, the degree to which each of these 
mechanisms contributes to the gene content evolution of 
X chromosomes remains unclear. 

Analyses of the evolutionary dynamics of sex-biased genes 
in Drosophila have shown that sex-biased genes evolve more 
rapidly than unbiased genes, which is expected because sexu- 
ally dimorphic traits contribute to reproductive success 
(Ellegren and Parsch 2007). Interestingly, this trend is driven 
by the evolution of male-biased genes, which have higher 
turnover (birth/death) rates and evolve faster at the sequence 
and expression levels than female-biased genes (Meikeljohn 
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et al. 2003; Zhang and Parsch 2005; Proschel et al. 2006; 
Sawyer et al. 2007; Zhang et al. 2007). Though some evidence 
shows that these differences are due to stronger positive 
selection acting on male-biased genes (Zhang and Parsch 
2005; Proschel et al. 2006; Sawyer et al. 2007), another 
factor may be that female-biased genes are evolutionarily con- 
strained due to greater functional pleiotropy (Parisi et al. 2004; 
Zhang et al. 2007; Mank et al. 2008; Mank and Ellegren 2009; 
Meikeljohn and Presgraves 2012). 

A major limitation of previous studies is the lack of a clear 
definition of a sex-biased gene (Ellegren and Parsch 2007; 
Meisel 2011). In particular, studies have used a wide range 
of experimental designs, analyses, and thresholds to classify 
sex-biased genes in Drosophila. Nearly all analyses of sex- 
biased genes in Drosophila are based on in-house microarray 
experiments (Meikeljohn et al. 2003; Ranz et al. 2003; 
Mclntyre et al. 2006; Sturgill et al. 2007; Zhang et al. 2007; 
Jiang and Machado 2009), in which conditions for raising, 
collecting, and performing experiments on flies varied 
widely. Moreover, in the majority of cases, sex-biased genes 
were identified solely by comparing expression between 
whole males and whole females (Meikeljohn et al. 2003; 
Ranz et al. 2003; Mclntyre et al. 2006; Sturgill et al. 2007; 
Zhang et al. 2007; Jiang and Machado 2009), whereas in a 
few, comparisons between ovary and testis tissue expression 
were also used (Parisi et al. 2003; Proschel et al. 2006). 
To identify sex-biased genes from these comparisons, some 
investigators used a 2-fold cutoff for male/female expression 
ratios (Parisi et al. 2003; Proschel et al. 2006), whereas others 
used a variety of statistical approaches with different assump- 
tions and significance thresholds (Meikeljohn et al. 2003; Ranz 
et al. 2003; Mclntyre et al. 2006; Sturgill et al. 2007; Jiang 
and Machado 2009). Because detection of sex-biased gene 
expression can be affected by differences in experimental 
conditions (Wyman et al. 2010) and statistical analyses, as 
well as by variation in biological factors, such as genotype 
(Jin et al. 2001; Gibson et al. 2004) and age (Jin et al. 
2001), it is not surprising that estimates of the fraction of 
genes with sex-biased expression in Drosophila range from 
10% to 91%. 

We investigated sex-biased transcriptome evolution be- 
tween D. melanogaster and D. pseudoobscura, two 
Drosophila species that diverged approximately 21^6 Ma 
(Beckenbach et al. 1993). In particular, the acquisition of an 
additional sex chromosome, XR, in D. pseudoobscura approxi- 
mately 8-12 Ma (Richards et al. 2005) enabled us to contrast 
gene content evolution on sex chromosomes of different 
ages. Because classification of sex-biased genes is fundamen- 
tal to studying their evolution, we first compared several 
independent metrics for identifying male- and female-biased 
genes. Then, limiting our analysis to genes that were classified 
as sex biased by multiple metrics, we studied the genomic 
properties and evolutionary dynamics of sex-biased transcrip- 
tome expression in Drosophila. 



Materials and Methods 

Gene Expression Data 

Paired-end RNA-seq reads from D. melanogaster whole male, 
whole female, carcass (mixed males and females), male head, 
female head, testis, accessory gland, and ovary tissues were 
downloaded from the modENCODE site at http://www.mod 
encode.org/. Paired-end RNA-seq reads for whole male, 
whole female, male carcass, female carcass, testis, accessory 
gland, and ovary tissues in D. pseudoobscura were obtained 
from Kaiser et al. (201 1), and comparable data for male and 
female head tissues in D. pseudoobscura were downloaded 
from NCBI's sequence read archive (accession numbers 
SRX016182 and SRX016183). 

Bowtie2 (Langmead et al. 2009) was used to align 
reads to transcript sequences of D. melanogaster and 
D. pseudoobscura, using annotation files downloaded 
from http://www.flybase.org (version 3 for both species). 
Transcript abundances were calculated using eXpress 
(Roberts and Pachter forthcoming), which outputs read 
counts and the number of fragments per kilobase of exon 
per million fragments mapped (FPKM; Trapnell et al. 2010). 
Quantile normalization of FPKMs was performed using the 
affy package of Bioconductor in the R software environment 
(R Development Core Team 2009). By comparing distributions 
of quantile-normalized FPKMs for exons to those for introns 
and intergenic sequences (obtained using the procedure out- 
lined for exons), we established cutoffs for expressed tran- 
scripts of FPKM = 1 and FPKM = 4 in D. melanogaster and 
D. pseudoobscura, respectively. Though different, these 
cutoffs resulted in similar proportions of expressed genes per 
tissue in the two species. To enable comparison between 
the two species, we scaled D. pseudoobscura FPKMs to 
those of D. melanogaster. 

Identification of Sex-Biased Genes 

To identify sex-biased genes in D. melanogaster and 
D. pseudoobscura, we applied several independent metrics. 
First, we simply compared the absolute expression of genes 
between whole males and whole females, and between testes 
to ovaries, and obtained genes for which the FPKM was 2-fold 
higher in one tissue than in the other. In addition to this naive 
approach, we also identified genes differentially expressed be- 
tween whole males and whole females, as well as between 
testes and ovaries, with three widely used tools: CuffDiff 
(Trapnell et al. 2010), DESeq (Anders 2010), and edgeR 
(Robinson et al. 2009). CuffDiff takes a nonparametric, 
annotation-guided approach to estimate the means and vari- 
ances of transcripts' FPKMs in different conditions, using 
Student's Mests to identify differentially expressed transcripts 
(Trapnell et al. 2010). In contrast, DESeq and edgeR estimate 
the means and variances of raw read counts under a negative 
binomial distribution and use exact tests to identify differen- 
tially expressed transcripts. The main difference between 
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DESeq and edgeR is that they use different statistical 
approaches to estimate variance (Robinson et al. 2009; 
Anders 2010). 

As a last metric, we obtained all sex tissue-specific genes 
(those primarily expressed in testes, accessory glands, or 
ovaries) by using the following formula to calculate the 
tissue specificity index, x, for each gene: 

1 _ log Ej 



where N is the number of tissues, E, is the expression in tissue 
/, and Emax is the maximum expression of the gene in all tissues 
(Yanai et al. 2005; Larracuente et al. 2008). x ranges from 0 
to 1, with larger x values indicating greater tissue specificity. 
We applied a cutoff of r > 0.9 to obtain highly tissue-specific 
genes (Larracuente et al. 2008). Because genes that are 
expressed in multiple tissues are more likely to have complex 
functions than those expressed in a single tissue (McShea 
2000), X can also be used to approximate the pleiotropy of 
genes (Mank et al. 2008; Mank and Ellegren 2009; Meisel 
2011; Meikeljohn and Presgraves 2012). Therefore, in our 
analyses, we assumed that broadly expressed genes (those 
with low X values) are more pleiotropic than narrowly 
expressed genes (those with high x values). 

Evolution of Sex-Biased Genes 

We obtained orthologs for sex-biased and unbiased genes 
from Drosophila ortholog tables downloaded from http:// 
www.flybase.org, which we supplemented with genes from 
reciprocal best Basic Local Alignment Search Tool (BLAST, 
Altschul et al. 1990) and BLAST-Like Alignment Tool (BLAT, 
Kent 2002) searches. Our turnover analysis was conservative 
in that, even if a one-to-one ortholog could not be identified 
for a gene in the sister species, it was not considered to be 
absent unless it had no potential orthologs in the genome of 
that species. Our analysis of expression divergence was also 
conservative; we only considered a gene as having gained/lost 
its sex-biased expression if it was identified as sex biased in one 
species by multiple metrics (n = 2 or n = 4; see main text for 
details) and not identified as sex biased by any metrics in the 
sister species. 

We downloaded D. melanogaster and D. pseudoobscura 
CDS sequences from http://www.flybase.org and aligned 
orthologous genes with MACSE (Ranwez et al. 2011), 
which performs in-frame alignments of protein-coding 
sequences. PAML (Yang 2007) was used to estimate the pair- 
wise KJKs from these alignments. 

Statistical Analyses 

We used Wilcoxon tests to compare the overall expression 
level, sex-biased expression, and tissue specificity among 
chromosomes, as well as to compare the KJKs of unbiased. 



male-biased, and female-biased genes. The expected number 
of male-/female-biased genes on a particular chromosome 
was determined by multiplying the total number of male-/ 
female-biased genes in the genome by the proportion of all 
genes on that chromosome, tests were used to determine 
whether the observed and expected numbers of male- and 
female-biased genes on sex chromosomes were significantly 
different. All statistical analyses were performed in the R soft- 
ware environment (R Development Core Team 2009). 

Results 

identification of Sex-Biased Genes in D. melanogaster 
and D. pseudoobscura 

We applied nine independent metrics to identify genes 
with sex-biased expression in D. melanogaster and 
D. pseudoobscura. First, we obtained genes differentially ex- 
pressed in whole males and whole females, as well as in testes 
and ovaries, by four analyses: comparison of absolute expres- 
sion levels between male and female tissues (2-fold cutoff), 
CuffDiff (Trapnell et al. 2010), DESeq (Anders 2010), and 
edgeR (Robinson et al. 2009). Additionally, we ascertained 
genes expressed primarily in the sex tissues of males (testis 
and accessory gland) or females (ovary) by calculating breadths 
of gene expression (see Materials and Methods for details). 

Depending on the metric used, 6.9-28.1% of genes 
in D. melanogaster and 6.6-24.1% of genes in 
D. pseudoobscura are classified as male biased, and 
0.3-28.5% of genes in D. melanogaster and 0.2-35.9% of 
genes in D. pseudoobscura are classified as female biased 
(table 1; see supplementary tables S1-S4, Supplementary 
Material online, for lists of genes and metrics). In total, ap- 
proximately 42.8% of genes in D. melanogaster and 31 % of 
genes in D. pseudoobscura are classified as male biased by at 
least one metric, and 34.8% of genes in D. melanogaster and 
46% of genes in D. pseudoobscura are classified as female 
biased by at least one metric. However, the overlap of genes 
satisfying multiple definitions is small, particularly for fe- 
male-biased genes (fig. 1 and supplementary tables S1-S4, 
Supplementary Material online). There are 339 and 493 
genes identified as male biased by all metrics in D. melanoga- 
ster and D. pseudoobscura, respectively, only 16 genes identi- 
fied as female biased by all metrics in D. melanogaster, and 
not a single gene identified as female biased by all metrics in 
D. pseudoobscura. Thus, male-biased genes are identified 
more robustly than female-biased genes. 

In total, more than 75% of genes in each species were 
classified as sex biased by at least one metric. To determine 
whether sets of genes identified by different metrics have 
unique properties, we compared the proportions of X-linked 
male- and female-biased genes among metrics. This analysis 
revealed that genes identified by different metrics tend 
to have different genomic distributions (supplementary 
table S5, Supplementary Material online). Hence, rather 
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Table 1 

Sex-Biased Genes Identified by Various Metrics 



Male-Biased Genes Female-Biased Genes 





D. melanogaster 


D. pseudoobscura 


D. melanogaster 


D. pseudoobscura 


Body, 2-fold 


3,871 


3,148 


3,512 


5,737 


Body, CuffDiff 


1,081 


1,980 


225 


188 


Body, DESeq 


948 


1,054 


39 


24 


Body, edgeR 


1,635 


1,964 


89 


55 


Gonads, 2-fold 


3,864 


3,851 


3,927 


4,826 


Gonads, CuffDiff 


3,713 


2,201 


1,559 


295 


Gonads, DESeq 


1,690 


1,348 


37 


30 


Gonads, edgeR 


2,450 


2,664 


105 


96 


Sex tissue-specific 


2,220 


2,412 


413 


408 




Fig. 1. — Comparison of four metrics for classifying sex-biased genes. Venn diagrams of numbers of male-biased (left) and female-biased 
(right) genes identified by four comparisons of whole male and whole female tissues in D. melanogaster (top) and D. pseudoobscura (bottom). 



than choosing a specific nnetric(s), we defined sex-biased 
genes as those identified by any combination of at least 
n metrics. To investigate the relationship between n and 
degree of sex-biased gene expression, we plotted n against 
whole male/whole female expression (fig. 2). A key observa- 
tion is that variance in sex-biased expression increases with n, 
illustrating the strong positive correlation between the 
number of metrics used and the degree of sex-biased effect 
of genes identified. Moreover, as n increases, the distance 
between genes classified as male- and female-biased also in- 
creases. As expected, n for male-biased genes is positively 



correlated to whole male/whole female expression (p = 0.80 
for D. melanogaster and p = 0.83 for D. pseudoobscura), 
whereas n for female-biased genes is negatively correlated 
to whole male/whole female expression (p = -0.60 for 
D. melanogaster and p = -0.22 for D. pseudoobscura). In 
general, the strengths of these correlations are consistent 
with the robustness of sex-biased gene identification (fig. 1 
and supplementary tables S1-S4, Supplementary Material 
online). 

In choosing n for defining sex-biased genes, we considered 
that requiring support from too few metrics would result in 
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nnany false positives (Type I errors), whereas requiring support 
fronn too nnany metrics would produce many false negatives 
(Type II errors). For this reason, we compared the numbers 
of overlapping and nonoverlapping sex-biased genes for dif- 
ferent n and selected two separate thresholds of n = 2 and 
n = 4. Genes that were classified as sex biased by at least 
any of two or four metrics were assigned to the n = 2 and 
n = 4 categories, respectively. Thus, n = 4 is a subset of n = 2. 
For n = 2, there are 4,042 male- and 3,038 female-biased 
genes in D. melanogaster, and 3,512 male- and 3,123 
female-biased genes in D. pseudoobscura. For n = 4, there 
are 2,281 male- and 331 female-biased genes in 
D. melanogaster, and 2,463 male- and 136 female-biased 
genes in D. pseudoobscura. 



D. ps&udoobscisra 



CL 

10 



D3 
O 

-5 



8 ^ 



U 



1 2 3 4 5 6 7 8 

n 



0 1 



r — r — r — r — r — r 
23456769 



Fig. 2. — Relationship between n and sex-biased gene expression. 
Plots show correlations between the number of metrics (n) that identify 
genes as male biased (blue) or female biased (red) and the whole male/ 
whole female expression ratios of those genes in D. melanogaster (left) and 
D. pseudoobscura (right). 



Chromosomal Distribution of Sex-Biased Genes in 
D. melanogaster and D. pseudoobscura 

Studies in Drosophila have revealed that there is an under- 
representation of male-biased genes (demasculinization) and 
over-representation of female-biased genes (feminization) 
on X chromosomes (Parisi et al. 2003; Sturgill et al. 2007). 
To investigate the extent of demasculinization and fem- 
inization in our data, we used four separate measures. First, 
we compared the absolute expression between genes on 
X chromosomes and autosomes in different tissues (fig. 3). 
Consistent with previous studies, we found that expression 
of genes in whole males, testes, and accessory glands is 
lower on the X chromosome in D. melanogaster and on 
both the XL and XR chromosomes in D. pseudoobscura, 
than on autosomes in either species, supporting demasculini- 
zation of X chromosomes. Additionally, the expression of 
genes in whole females and ovaries is significantly higher on 
the X chromosome than on autosomes in D. melanogaster, 
supporting feminization of the D. melanogaster X, though 
there is no evidence for feminization of either XL or XR in 
D. pseudoobscura. 

Second, we compared the contribution of sex-specific 
tissue to whole body expression (testis/whole male, accessory 
gland/whole male, and ovary/whole female expression) 
on X chromosomes and autosomes in D. melanogaster 
and D. pseudoobscura (fig. 4). In both species, we found 
evidence for demasculinization of X chromosomes when 
comparing testis/whole male contributions and feminiza- 
tion when comparing ovary/whole female contributions. 
Surprisingly, however, accessory gland tissue contributions 




-5 0 5 10 -5 0 5 10 -5 0 5 10 -a 0 5 10 -5 0 S 10 
log{FPKM} log(FR<M) lc3§[FPKM) log(FPKM) 1og{FPKM) 



Fig. 3.— Absolute expression on X chromosomes and autosomes in different tissues. Density plots of absolute expression for X chromosomes 
and autosomes in different tissue of D. melanogaster (top) and D. pseudoobscura (bottom). Each panel corresponds to a different tissue, with autosomal 
densities depicted in blue, ancestral X chromosome (D. melanogaster X and D. pseudoobscura XL) densities depicted in red, and D. pseudoobscura 
XR densities depicted in orange. Medians are given in the top left corner of each panel, with asterisks corresponding to P< 0.05 (*), P< 0.01 (**), and 
P< 0.001 (***). 



Genome Biol. EvoL 4(1 1):1 189-1200. doi:10.1093/gbe/evs093 Advance Access publication October 23, 2012 



1193 



Assis et al. 



GBE 




Fig. 4. — Contributions of sex-specific tissue expression on X chromosomes and autosomes. Density plots of sex-specific tissue/whole body expression in 
D. melanogaster (top) and D. pseudoobscura (bottom). Each panel corresponds to a different sex-specific tissue, with autosomal densities depicted in blue, 
ancestral X chromosome (D. melanogaster X and D. pseudoobscura XL) densities depicted in red, and D. pseudoobscura XR densities depicted in orange. 
Medians are given in the top left corner of each panel, with asterisks corresponding to P< 0.05 (*), P< 0.01 (**), and P< 0.001 (***). 



were significantly higher on the X chromosome in D. melano- 
gaster and on both XL and XR in D. pseudoobscura than on 
autosomes in either species, providing evidence for masculin- 
ization of Drosophila X chromosomes. 

Third, we compared sex-biased expression of whole 
body (whole male/whole female) and gonad (testis/ovary) 
tissues on X chromosomes and autosomes (fig. 5). 
Consistent with previous studies, we found that whole 
male/whole female and testis/ovary expression ratios are sig- 
nificantly higher on autosomes than on X chromosomes in 
D. melanogaster and D. pseudoobscura, providing evidence 
for demasculinization and/or feminization of the X chromo- 
somes in both species. 

As a final approach, we compared observed and 
expected numbers of male- and female-biased genes on 
D. melanogaster and D. pseudoobscura X chromosomes 
(fig. 6). We calculated expected numbers by assuming 
that sex-biased genes are randomly distributed across 
chromosomes and that their representation on a particular 
chromosome is proportional to the actual number of genes 
on that chromosome. For n = 2, the observed numbers 
of male-biased genes on D. melanogaster X and on 
D. pseudoobscura XL and XR are significantly lower than 
expected, and the observed numbers of female-biased 
genes on D. melanogaster X and on D. pseudoobscura XL 
and XR are significantly higher than expected, supporting 
demasculinization and feminization of all Drosophila 
X chromosomes. However, for n = 4, demasculinization is 
not observed for D. pseudoobscura XL, and the significance 




Sex-specHic 



-0.31 

J 








0.26 
-0.76*** 

-0,67*** £K 









-6^-2 0246 -6 ^^2 0246 
log(male / female FPKM) log(testis / ovary FPKM) 

Fig. 5. — Sex-biased expression on X chromosomes and autosomes. 
Density plots of whole male/whole female (left) and testis/ovary (right) 
expression ratios in D. melanogaster (top) and D. pseudoobscura 
(bottom). Autosomal densities are depicted in blue, ancestral X chromo- 
some (D. melanogaster X and D. pseudoobscura XL) densities are depicted 
in red, and D. pseudoobscura XR densities are depicted in orange. Vertical 
dashed lines indicate equal male and female expression. Medians are given 
in the top left corner of each panel, with asterisks corresponding to 
P< 0.05 (*), P< 0.01 (**), and P< 0.001 (***). 
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Fig. 6. — Representation of sex-biased genes on Drosophila X chromosomes. Proportions of observed/expected numbers of male-biased (blue) and 
female-biased (red) genes on D. melanogasterX, D. pseudoobscura XL, and D. pseudoobscura XR. The number of metrics (n) used in the classification of 
sex-biased genes for each analysis is indicated on the x axis. Horizontal dashed lines indicate equal observed and expected numbers of genes, and asterisks 
indicate P< 0.05 (*), P< 0.01 (**), and P< 0.001 (***). 



of feminization of D. pseudoobscura XR is weaker than 
for n = 2. 

Though most approaches revealed a general trend of 
demasculinization and feminization of Drosophila X chromo- 
somes, aside from more expressed X-linked genes in 
D. melanogaster ovaries, the total numbers of expressed 
genes in different tissues do not differ between the X chromo- 
somes and autosomes. Therefore, demasculinization does 
not simply result from a lower proportion of X-linked genes 
expressed in testis, as might be expected if X chromosomes 
are transcriptionally silenced during spermatogenesis or if 
sexually antagonistic selection removes individual male-biased 
genes. Rather, lower male-biased expression appears to be 
a property of most genes residing on X chromosomes. 

Evolution of Sex-Biased Genes between D. melanogaster 
and D. pseudoobscura 

To study the evolution of sex-biased genes, we obtained 
pairs of orthologous genes in D. melanogaster and 
D. pseudoobscura and calculated the turnover, expression di- 
vergence, and sequence evolutionary rates of genes identified 
as male and female biased in the two species. Turnover of a 
sex-biased gene was inferred when that gene did not have an 
ortholog in the sister species (see Materials and Methods for 
details), and the rate of turnover was estimated by dividing the 
number of such genes by the total number of sex-biased 
genes. Comparison of sex-biased gene turnover rates con- 
firmed previous findings that, in D. melanogaster and 
D. pseudoobscura, male-biased genes have much higher turn- 
over rates than female-biased genes (for n = 2 and n = 4). 
Although nearly all female-biased genes are present in both 
species, approximately 19.6-37.6% of male-biased genes in 
either D. melanogaster or D. pseudoobscura lack orthologs in 
the other species (fig. 7 A). 

Expression divergence was inferred when a gene was clas- 
sified as male or female biased (n = 2 or n = 4) in one species 
and not classified as male/female biased by any metrics in the 
sister species. The rate of evolution by expression divergence 
was estimated by dividing the number of genes with 



expression divergence by the total number of sex-biased 
genes containing orthologs in the sister species. Comparison 
of the rates between male- and female-biased genes revealed 
opposite patterns in D. melanogaster and D. pseudoobscura 
(fig. 7B). 

In D. melanogaster, higher proportions of male- than 
female-biased genes do not have the same sex bias in 
D. pseudoobscura. On the other hand, in D. pseudoobscura, 
more female- than male-biased genes do not have the same 
sex bias in D. melanogaster. Interestingly, this difference 
seems to be due to contrasting rates of switches between 
male- and female-biased states in the two species (table 2). 
In particular, from D. melanogaster to D. pseudoobscura, 
male-to-female-biased switches occur 4.2-5.6 times faster 
than female-to-male-biased switches. In contrast, from 
D. pseudoobscura to D. melanogaster, female-to-male- 
biased switches occur 2.7-4.7 times faster than male-to- 
female-biased switches. Thus, from D. melanogaster to 
D. pseudoobscura, changes in sex-biased states tend to be 
in the male-to-female direction, whereas the opposite pattern 
is observed from D. pseudoobscura to D. melanogaster 

To study the sequence divergence of sex-biased genes, we 
calculated sequence evolutionary rates (KJKs) of unbiased, 
male-biased, and female-biased genes in D. melanogaster 
and D. pseudoobscura. In both species, for n = 2 and n = 4, 
sex-biased genes evolve much faster than unbiased genes at 
the DNA sequence level (P< 0.001). Comparison of male- 
and female-biased genes revealed that male-biased genes 
evolve significantly faster than female-biased genes (fig. 7C, 
P< 0.001), which evolve at similar rates as unbiased genes. 
Thus, we confirmed previous findings that faster evolution of 
sex-biased genes is primarily due to the rapid evolution of 
male-biased genes. 

Faster evolution of male-biased genes is often attributed 
to stronger positive selection acting on male- than on fe- 
male-biased genes (Zhang and Parsch 2005; Proschel et al. 
2006; Sawyer et al. 2007). However, another possible factor 
is that female-biased genes are more pleiotropic, restricting 
their adaptive potential (Parisi et al. 2004; Zhang et al. 2007). 
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Fig. 7. — Evolution of sex-biased genes in Drosophila. Rates of D. melanogaster (top) and D. pseudoobscura (bottom) male-biased (blue) and fe- 
male-biased (red) gene tumover (A), expression divergence (B), and sequence divergence (0- The number of metrics (n) used in the classification of 
sex-biased genes for each analysis is indicated on the x axis. 



Table 2 

Rates of Switches between Male- and Female-Biased Gene States 





n = 2 






n = 4 










F-^M 


D. melanogaster 


0.187 


0.045 


0.067 


0.012 


D. pseudoobscura 


0.051 


0.139 


0.014 


0.066 



This hypothesis is supported by the observation that most 
male-biased expression in D. melanogaster results from 
genes expressed in a sex-specific tissue(s), whereas 
female-biased genes tend to be broadly expressed (Parisi 
et al. 2004; Sturgill et al. 2007; Meikeljohn and Presgraves 
2012). Examination of the contribution of individual tissues 
to sex-biased expression revealed that only testis and acces- 
sory gland expression levels are positively correlated to 
sex-biased expression (P<0.05), illustrating the influence of 
male sex tissue-specific expression to overall male-biased gene 
expression (supplementary fig. S1, Supplementary Material 
online). To further investigate the pleiotropy hypothesis, we 
compared tissue specificity (x) to whole male/whole female 
expression ratios in D. melanogaster and D. pseudoobscura 
(fig. 8). This analysis shows that, in both species, female-biased 
genes tend to have smaller x values than male-biased genes, 
supporting the hypothesis that female-biased genes exhibit 
greater pleiotropy. Thus, female-biased genes are expected 
to be under greater selective constraint than male-biased 
genes, likely contributing to their slower evolution. 




a2 OA 0.6 0.8 1.0 

t 



0.2 0.4 0.6 0.8 1.0 



Fig. 8. — ^Tissue specificity of sex-biased genes. Correlations between 
whole male/whole female expression ratios and x in D. melanogaster {left) 
and D. pseudoobscura (right). Unbiased genes are represented by gray 
dots, male-biased genes by blue dots, and female-biased genes by red 
dots {n = 2). Horizontal dashed lines indicate equal male and female 
expression. 



To determine whether gene turnover or switches in 
sex-biased expression contribute to demasculinization of 
X chromosomes, we compared the turnover and expression 
evolutionary rates of different Muller elements, which repre- 
sent orthologous chromosomal segments, in D. melanogaster 
and D. pseudoobscura (fig. 9). We found that, in each species, 
male-biased genes on Muller-A, the ancestral X chromosome, 
are least likely to have orthologs in the other species 
(P<0.01). Furthermore, genes that are male biased in 
D. melanogaster are most likely to show female-biased expres- 
sion in D. pseudoobscura, and genes that are female biased 
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Fig. 9. — ^Turnover and expression evolution of sex-biased genes among chromosomes. Proportions of male- and female-biased genes on each Muller 
element of D. melanogaster (top) and D. pseudoobscura (bottom) that, in the other species, have conserved expression, are unbiased, have opposite 
sex-biased expression, or do not have orthologs. 



in D. pseudoobscura are most likely to display male-biased 
expression in D. melanogaster, when they are located 
on Muller-D, which is autosomal chromosome 3L in 
D. melanogaster and chromosome XR in D. pseudoobscura. 
Thus, both turnover in gene content and switches in expres- 
sion bias contribute to the deficiency of male-biased genes on 
Drosophila X chromosomes. 

Discussion 

Classification of Sex-Biased Genes Is Complex 

The continuous nature of sex-biased gene expression makes 
it difficult to consistently identify genes with the greatest 
sex-biased effects, complicating meaningful comparisons 
among taxa and hindering our understanding of sex-biased 
transcriptome evolution. In this study, we compared nine 
independent metrics for identifying sex-biased genes in 
D. melanogaster and D. pseudoobscura. We found that dif- 
ferent metrics yielded vastly different numbers and sets of 
sex-biased genes in both species, often with little overlap in 
the genes identified across metrics. However, this is not too 
surprising, given that the metrics are based on different bio- 
logical and/or statistical models and have different assump- 
tions, statistical methodologies, and false positive/negative 
rates. Because analysis of the sex-biased genes identified by 
different metrics did not reveal an optimal metric or set of 
metrics, we chose to define a sex-biased gene based on its 
support by a minimum number of metrics, n. We compared 
numbers of sex-biased genes for different n, and, to minimize 



Type I and Type II errors, we chose two separate cutoffs for the 
number of metrics supporting each gene: n = 2 and n = 4. 
Because male-biased genes are more robustly identified 
across metrics, each of these cutoffs yielded more male- 
than female-biased genes, with n = 2 resulting in over ten 
times more genes classified as female biased than n = 4. 
There are two possible reasons for this effect. First, there 
may be few female-biased genes, and, therefore, higher strin- 
gency (n = 4) enhances the identification of female-biased 
genes by eliminating Type I errors. On the other hand, fe- 
male-biased genes may be more difficult to identify, possibly 
because they are pleiotropic or expressed at lower levels than 
male-biased genes, leading to a general lack of female-biased 
genes among metrics. In this case, Type II errors are more 
problematic, perhaps making n = 2 a better cutoff for fe- 
male-biased genes. Thus, there was a twofold benefit to 
using both cutoffs in our analysis. It enabled us to compare 
conclusions reached by using different cutoffs, as well as to 
uncover well-supported patterns in our data. 

Drosophila X Chromosomes Are Demasculinized and 
Feminized 

Analysis of chromosomal gene content in Drosophila by 
four separate approaches revealed a general deficiency of 
male-biased genes and enrichment of female-biased genes 
on Drosophila X chromosomes. The exception to this result 
is accessory gland/whole male expression, which instead sup- 
ports masculinization of X chromosomes (fig. 4). Interestingly, 
demasculinization of X chromosomes is observed from overall 
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accessory gland expression patterns (fig. 3). To investigate this 
effect, we compared accessory gland/male carcass tissue 
expression on D. pseudoobscura X chromosomes and auto- 
somes. Unfortunately, because of the lack of male carcass ex- 
pression data for D. melanogaster, we were unable to perform 
this analysis in both species. In D. pseudoobscura, accessory 
gland/male carcass tissue expression was higher on autosomes 
than on X chromosomes, providing support for demasculiniza- 
tion of X chromosomes. Because male carcass tissue includes all 
individual tissues found in whole males, with the exception of 
sex-specific (testis and accessory gland) tissues, this result sug- 
gests that the masculinization pattern observed from accessory 
gland/whole male comparisons is driven by lower testis expres- 
sion on X chromosomes relative to autosomes. Thus, this pat- 
tern does not appear to be biologically interesting, and the 
general finding is that X chromosomes are demasculinized. 

Though D. pseudoobscura XR only became a sex chromo- 
some 8-12 Ma (Richards et al. 2005), it already shows evi- 
dence of demasculinization and feminization, revealing that 
gene content evolution occurs rapidly after acquisition of 
sex chromosomes. Interestingly, in D. pseudoobscura, demas- 
culinization is stronger on XR than on XL, the ancestral 
X chromosome, and feminization is stronger on XL than 
on XR. Because XR is a newer sex chromosome, this effect 
may be caused by demasculinization proceeding at a faster 
rate than feminization, presumably because male-biased 
genes evolve faster than female-biased genes. 

MSG! (Lifschytz and Lindsley 1972; Hense et al. 2007; 
Vibranovski et al. 2009), dosage compensation (Bachtrog 
et al. 2010), and sexually antagonistic selection (Rice 1984; 
Ellegren and Parsch 2003; Gurbich and Bachtrog 2008) 
have each been proposed to explain demasculinization of 
X chromosomes. Under the MSG! hypothesis, we would 
expect to observe demasculinization only in testis and 
whole male tissues, and we would also expect there to be 
fewer X-linked genes expressed in testes. Although we do 
observe demasculinization in testes and whole males, we 
also observe demasculinization in accessory glands when 
we consider absolute expression of genes. Moreover, 
there is no significant difference in the proportions of 
testis-expressed genes on X chromosomes and autosomes. 
These observations, together with the recent finding that 
MSG! may not operate in Drosophila (Meikeljohn et al. 
2011), suggest that MSG! is not responsible for demasculi- 
nization of Drosophila X chromosomes. In contrast, dosage 
compensation could cause demasculinization patterns in 
testes and whole males, as well as in accessory glands, for 
two separate reasons. Recent data suggest that testes 
in Drosophila completely lack dosage compensation 
(Meikeljohn et al. 2011), and generally lower expression 
levels of X-linked genes in testes are consistent with that 
finding (Meikeljohn and Presgraves 2012). In accessory 
glands, the dosage compensation machinery could interfere 
with the evolution of male-biased gene expression (Parisi 



et al. 2004). Additionally, dosage compensation would not 
reduce the proportion of genes expressed in either testes or 
accessory glands on the X chromosome, which is consistent 
with our findings. Sexually antagonistic selection, on the 
other hand, operates on a gene-by-gene basis. Thus, we 
would expect to observe fewer testis-expressed and more 
ovary-expressed genes on X chromosomes, which is not 
supported by our data. Instead, lower expression appears 
to be a general property of genes residing on Drosophila 
X chromosomes. Thus, our data support dosage compensa- 
tion as a major mechanism driving demasculinization of 
Drosophila X chromosomes. However, it is important to 
note that many independent forces likely contribute to sex 
chromosome evolution, and disentangling these forces may 
be impossible from a chromosome-wide perspective. 

Male-Biased Genes Evolve Faster than 
Female-Biased Genes 

Gonsistent with previous studies (Meikeljohn et al. 2003; 
Zhang and Parsch 2005; Proschel et al. 2006; Sawyer et al. 
2007; Zhang et al. 2007), our analyses revealed that 
male-biased genes have higher turnover and sequence evolu- 
tionary rates than female-biased genes in D. melanogaster and 
D. pseudoobscura. We also found that evolution by expression 
divergence is faster for male-biased genes in D. melanogaster 
and for female-biased genes in D. pseudoobscura. Further 
analysis revealed that this contrast is likely driven 
by species-specific differences in rates of switches between 
male-to-female-biased and female-to-male-biased expression. 
In particular, male-biased genes in D. melanogaster a^e more 
likely to display female-biased expression in D. pseudoobscura, 
and female-biased genes in D. pseudoobscura are more likely 
to display male-biased expression in D. melanogaster. The 
highest switch rates in D. pseudoobscura occur among fe- 
male-biased genes on chromosome XR, which is a recently 
acquired sex chromosome in D. pseudoobscura and is ortho- 
logous to autosomal chromosome 3L in D. melanogaster. The 
most likely explanation for this pattern is that many genes 
residing on the ancestral autosomal chromosome switched 
from male- to female biased after it became a sex chromo- 
some in the D. pseudoobscura lineage. Thus, differences in 
switch rates may contribute to demasculinization and femin- 
ization of D. pseudoobscura chromosome XR. 

The faster evolution of male-biased genes is generally 
attributed to stronger positive selection acting on male-biased 
relative to female-biased genes (Zhang and Parsch 2005; 
Proschel et al. 2006; Sawyer et al. 2007), which is supported 
by the higher KJK^ rates of male-biased genes. However, 
another factor that may contribute to this pattern is functional 
pleiotropy of female-biased genes, which would limit their 
evolutionary potential (Parisi et al. 2004; Zhang et al. 2007). 
The similarity of KJK^ in female-biased and unbiased genes is 
consistent with this hypothesis. Moreover, comparison of 
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tissue specificities of nnale- and female-biased genes revealed 
that female-biased genes are generally more broadly ex- 
pressed than male-biased genes. Thus, our findings support 
the hypothesis that pleiotropy of female-biased genes con- 
strains their evolution, contributing to their slower evolution 
relative to male-biased genes. 

Supplementary Material 

Supplementary figure S1 and tables S1-S5 are available 
at Genome Biology and Evolution online (http:/AAAAAA/.gbe. 
oxfordjournals.org/). 

Acknowledgments 

The authors thank three anonymous reviewers for their 
helpful comments. This work was supported by a Packard 
Fellowship to D.B. and NIH grants R01GM076007 and 
R01GM093182. 

Literature Cited 

Altschul SF, Gish W, Miller W, Myers EW, Upman DJ. 1990. Basic local 

alignment search tool. J Mol Biol. 5:403^10. 
Anders S. 2010. Analyzing RNA-Seq data with the "DESeq" package. 

Mol Biol. 43:1-17. 

Bachtrog D., Toda NRT, Lockton S. 2010. Dosage compensation and 
demasculinization of X chromosomes in Drosophila. Curr Biol. 20: 
1476-1481. 

Beckenbach AT, Wei YW, Liu H. 1993. Relationships in the Drosophila 
obscura species group, inferred from mitochondrial cytochrome 
oxidase II sequences. Mol Biol Evol. 10:619-634. 

Connallon J, Knowles L. 2005. Intergenomic conflict revealed by patterns 
of sex-biased gene expression. Trends Genet. 21 :495-499. 

Ellegren H, Parsch J. 2007. The evolution of sex-biased genes and 
sex-biased gene expression. Nat Rev Genet. 8:689-698. 

Gibson G, et al. 2004. Extensive sex-specific nonadditivity of gene expres- 
sion in Drosophila melanogaster. Genetics 167:1791-1799. 

Gurbich TA, Bachtrog D. 2008. Gene content evolution on the X chromo- 
some. Curr Opin Genet Dev. 18:493^98. 

Hense W, Baines JF, Parsch J. 2007. X chromosome inactivation during 
Drosophila spermatogenesis. PLoS Biol. 5:e273. 

Jiang Z, Machado CA. 2009. Evolution of sex-dependent gene expression 
in three recently diverged species of Drosophila. Genetics 183: 
1175-1185. 

Jin W, et al. 2001. The contributions of sex, genotype, and age to 
transcriptional variance in Drosophila melanogaster. Nature 29: 
389-395. 

Kaiser VB, Zhou Q, Bachtrog D. 2011. Non-random gene loss from the 
Drosophila miranda neo-Y chromosome. Genome Biol Evol. 3: 
1329-1337. 

Kent WJ. 2002. BLAT— the BLAST-like alignment tool. Genome Res. 12: 
656-664. 

Langmead B, Trapnell C, Pop M, Salzberg SL. 2009. Ultrafast and 

memory-efficient alignment of short DNA sequences to the human 

genome. Genome Biol. 10:R25. 
Larracuente AM, et al. 2008. Evolution of protein-coding genes in 

Drosophila. Trends Genet. 24:1 14-123. 
Lifschytz E, Undsley DL. 1972. The role of X-chromosome inactivation 

during spermatogenesis. Proc Natl Acad Sci USA. 69:182-186. 
Mank JE, Ellegren H. 2009. Are sex-biased genes more dispensable? 

Biol Lett. 5:409^12. 



Mank JE, Hultin-Rosenberg L, Zwahlen M, Ellegren H. 2008. Pleiotropic 

constraint hampers the resolution of sexual antagonism in vertebrate 

gene expression. Am Nat. 171:35-43. 
Mclntyre LM, et al. 2006. Sex-specific expression of alternative transcripts 

in Drosophila. Genome Biol. 7:R79. 
McShea DW. 2000. Functional complexity in organisms: parts as proxies. 

Biol Philos. 15:641-668. 
Meikeljohn CD, Landeen EL, Cook JM, Kingan SB, Presgraves DC. 201 1. 

Sex chromosome-specific regulation in the Drosophila male germline 

but little evidence for chromosomal dosage compensation or meiotic 

inactivation. PLoS Biol. 9:e1001126. 
Meikeljohn CD, Parsch J, Ranz JM, HartI DL. 2003. Rapid evolution 

of male-biased gene expression in Drosophila. Proc Natl Acad Sci 

USA. 17:9894-9899. 
Meikeljohn CD, Presgraves DC. 2012. Little evidence for demasculinization 

of the Drosophila X chromosome among genes expressed in the male 

germline. Genome Biol Evol. 4:895-904. 
Meisel RP. 201 1 . Towards a more nuanced understanding of the relation- 
ship between sex-biased gene expression and rates of protein-coding 

sequence evolution. Mol Biol Evol. 28:1893-1900. 
Parisi M, et al. 2003. Paucity of genes on the Drosophila X chromosome 

showing male-biased expression. Science 299:697-700. 
Parisi M, et al. 2004. A survey of ovary-, testis-, and soma-biased 

gene expression in Drosophila melanogaster adults. Genome Biol. 

5:R40. 

Proschel M, Zhang Z, Parsch J. 2006. Widespread adaptive evolution 
of Drosophila genes with sex-biased expression. Genetics 174: 
893-900. 

R Development Core Team. 2009. R: a language and environment for 
statistical computing. Vienna (Austria): R Foundation for Statistical 
Computing. 

Ranwez V, Harispe S, Delsuc F, Douzery EJP. 201 1 . MACSE: multiple align- 
ment of coding sequences accounting for frameshifts and stop 
codons. PLoS One 6(9):e22594. 

Ranz JM, Castillo-Davis CI, Meikeljohn CD, HartI DL. 2003. Sex-dependent 
gene expression and evolution of the Drosophila transcriptome. 
Science 300:1742^745. 

Rice WR. 1984. Sex chromosomes and the evolution of sexual dimorph- 
ism. Evolution 38:735-742. 

Richards S, et al. 2005. Comparative genome sequencing of Drosophila 
pseudoobscura: chromosomal, gene, and cis-element evolution. 
Genome Res. 15:1-18. 

Rinn JL, Snyder M. 2005. Sexual dimorphism in mammalian gene expres- 
sion. Trends Genet. 21:298-305. 

Roberts A, Pachter L. Forthcoming. Streaming fragment assignment for 
real-time analysis of sequencing experiments. Nat Methods. 

Robinson MD, McCarthy DJ, Smyth GK. 2009. edgeR: a bioconductor 
package for differential expression analysis of digital gene expression 
data. Bioinformatics 26:139-140. 

Sawyer SA, Parsch J, Zhang Z, HartI DL. 2007. Prevalence of positive 
selection among nearly neutral amino add replacements in 
Drosophila. Proc Natl Acad Sci USA. 104:6504-6510. 

Sturgill D, Zhang Y, Parisi M, Oliver B. 2007. Demasculinization of 
X chromosomes in the Drosophila genus. Nature 450:238-241. 

Trapnell C, et al. 2010. Transcript assembly and quantification by RNA-Seq 
reveals unannotated transcripts and isoform switching during cell 
differentiation. Nat Biotechnol. 28:511-515. 

Vibranovski MD, Lopes HF, Karr TL, Long M. 2009. Stage-specific expres- 
sion profiling of Drosophila spermatogenesis suggests that meiotic sex 
chromosome inactivation drives genomic relocation of testis-expressed 
genes. PLoS Genet. 5:e1 000731. 

Vicoso B, Charlesworth B. 2009. The deficit of male-biased genes on the 
D. melanogaster X is expression-dependent: a consequence of dosage 
compensation? J Mol Evol. 68:576-583. 



Genome Biol Evol 4(1 1):1 189-1200. doi:10.1093/gbe/evs093 Advance Access publication October 23, 2012 



1199 



Assis et al. 



GBE 



Wyman MJ, Agrawal AF, Rowe L. 2010. Condition-dependence of the 

sexually dimorphic transcriptome in Drosophila melanogaster. 

Evolution 64-66:1836-1848. 
Yanai I, et al. 2005. Genome-wide midrange transcription profiles reveal 

expression level relationships in human tissue specifications. 

Bioinformatics 21:650-659. 
Yang Z. 2007. PAML4: phylogenetic analysis by maximum likelihood. Mol 

Biol Evol. 24:1586-1591. 



Zhang Y, Sturgill D, Parisi M, Kumar S, Oliver B. 2007. Constraint and 

turnover in sex-biased gene expression in the genus Drosophila. 

Nature 450:233-237. 
Zhang Z, Parsch J. 2005. Positive correlation between evolutionary rate and 

recombination rate in Drosophila genes with male-biased expression. 

Mol Biol Evol. 22:1945-1947. 

Associate editor: Judith Mank 



1 200 Genome Biol. Evol. 4(1 1):1 189-1200. doi:10.1093/gbe/evs093 Advance Access publication October 23, 2012 



